##################################################################PLOTTING#########################################################
library(grid)
library(lattice)
library(latticeExtra)
trellis.par.set(strip.background=list(col="grey95"))
########################################################## 10 minute station annual ###############################################
# RAW data
#X/Y/Z
stationplotraw <- xyplot(stationplot[,c(1:8)], layout=c(1,8), main="Raw station data",
scales=list(y=list(cex=0.7, rot=0)),
col="black", lwd=0.8,
xlim=as.POSIXct(c("2014-06-01 00:00:00","2016-09-01 00:00:00")),
ylab=expression("mm","Wm"^-2,~degree~C,~degree~C,~degree~C,"ms"^-2,"ms"^-2,"ms"^-2),more=TRUE)
update(stationplotraw, panel = function(...) {
panel.abline(h = 0, lty = "dotted", col = "red")
panel.xyplot(...)
})
#X/Y/Z
stationplotraw <- xyplot(stationplot[,c(1:3)], layout=c(1,3), main="Raw station data",
scales=list(y=list(cex=0.7, rot=0)),
col="black",
xlim=as.POSIXct(c("2015-07-02 00:00:00","2015-07-31 00:00:00")),
ylab=expression("ms"^-1,"ms"^-1,"ms"^-1),more=TRUE)
update(stationplotraw, panel = function(...) {
panel.abline(h = 0, lty = "dotted", col = "red")
panel.xyplot(...)
})
a=xyplot((onejulyzoowindow[,c(1:3)]), layout=c(1,3), main="One second data against averaged 10 minute data over one week",
scales=list(y=list(cex=0.7, rot=0)),
key=list(columns=2,text=list(lab=c("ten minute average","one second")),lines=list(col=c("black","grey"))),
col="grey",
ylab=expression("ms"^-2,"ms"^-2,"ms"^-2),
xlab= "Days (2016)",
xlim=as.POSIXct(c("2016-07-03 00:00:00.000","2016-07-06 00:00:00.000")), more=TRUE)
bg <- trellis.par.get("background")
bg$col <- "transparent"
trellis.par.set("background", bg)
b=xyplot((stationplot[,c(1:3)]), layout=c(1,3),
scales=list(y=list(cex=0.7, rot=0)),
col="black",
ylab=expression("ms"^-2,"ms"^-2,"ms"^-2),
)
a + as.layer(b, x.same=TRUE,axes = TRUE)
#Weekly averages
#first convert to xts data so can use daily/monthly/yearly functions, which will be used to calculate means at different time periods
library(xts)
stationannual <- apply.weekly((as.xts(stationplot)), colMeans)
#X
plot.new()
a=xyplot((stationannual[,c(1:8)]), layout=c(1,8), main="Averaged daily weekly from mid-2014 to mid-2016",
scales=list(y=list(cex=0.7, rot=0)),
key=list(columns=2,text=list(lab=c("2014-2015","2015-2016")),lines=list(col=c("darkred","blue"))),
col="darkred",
ylab=expression("mm","Wm"^-2,~degree~C,~degree~C,~degree~C,"ms"^-1,"ms"^-1,"ms"^-1),
xlim=as.POSIXct(c("2014-08-01 00:00:00","2015-08-01 00:00:00")), more=TRUE)
bg <- trellis.par.get("background")
bg$col <- "transparent"
trellis.par.set("background", bg)
b=xyplot((stationannual[,c(1:8)]), layout=c(1,8),
scales=list(y=list(cex=0.7, rot=0)),
col="blue",
ylab=expression("mm","Wm"^-2,~degree~C,~degree~C,~degree~C,"ms"^-1,"ms"^-1,"ms"^-1),
xlim=as.POSIXct(c("2015-08-01 00:00:00","2016-08-01 00:00:00")))
a + as.layer(b, x.same=FALSE,axes = FALSE)
############################################################ 10 minute seasonal ###################################################
#daily averages
stationseasonal <- apply.daily((as.xts(stationplot)), colMeans)
# SPRING
plot.new()
a=xyplot((stationseasonal[,c(1:8)]), layout=c(1,8), main="Spring",
scales=list(y=list(cex=0.7, rot=0)),
key=list(columns=2,text=list(lab=c("2014-2015","2015-2016")),lines=list(col=c("darkred","blue"))),
col="darkred",
ylab=expression("mm","Wm"^-2,~degree~C,~degree~C,~degree~C,"ms"^-1,"ms"^-1,"ms"^-1),
ylim=list(c(0, 0.04),c(-0.08,0),c(-0.1,-0.02),c(10,25),c(10,25),c(10,25),c(0,400),c(0,0.4)),
xlim=as.POSIXct(c("2014-09-01 00:00:00","2014-12-01 00:00:00")), more=TRUE)
bg <- trellis.par.get("background")
bg$col <- "transparent"
trellis.par.set("background", bg)
b=xyplot((stationseasonal[,c(1:8)]), layout=c(1,8),
scales=list(y=list(cex=0.7, rot=0)),
col="blue",
ylab=expression("mm","Wm"^-2,~degree~C,~degree~C,~degree~C,"ms"^-1,"ms"^-1,"ms"^-1),
xlim=as.POSIXct(c("2015-09-01 00:00:00","2015-12-01 00:00:00")))
a + as.layer(b, x.same=FALSE,axes = FALSE)
# SUMMER
a=xyplot((stationseasonal[,c(1:8)]), layout=c(1,8), main="Summer",
scales=list(y=list(cex=0.7, rot=0)),
key=list(columns=2,text=list(lab=c("2014-2015","2015-2016")),lines=list(col=c("darkred","blue"))),
col="darkred",
ylab=expression("mm","Wm"^-2,~degree~C,~degree~C,~degree~C,"ms"^-1,"ms"^-1,"ms"^-1),
ylim=list(c(-0.02, 0.04),c(-0.04,0.01),c(-0.1,0.02),c(10,30),c(10,30),c(10,30),c(0,600),c(0,0.4)),
xlim=as.POSIXct(c("2014-12-01 00:00:00","2015-03-01 00:00:00")), more=TRUE)
bg <- trellis.par.get("background")
bg$col <- "transparent"
trellis.par.set("background", bg)
b=xyplot((stationseasonal[,c(1:8)]), layout=c(1,8),
scales=list(y=list(cex=0.7, rot=0)),
col="blue",
ylab=expression("mm","Wm"^-2,~degree~C,~degree~C,~degree~C,"ms"^-1,"ms"^-1,"ms"^-1),
xlim=as.POSIXct(c("2015-12-01 00:00:00","2016-03-01 00:00:00")))
a + as.layer(b, x.same=FALSE,axes = FALSE)
# AUTUMN
a=xyplot((stationseasonal[,c(1:8)]), layout=c(1,8), main="Autumn",
scales=list(y=list(cex=0.7, rot=0)),
key=list(columns=2,text=list(lab=c("2014-2015","2015-2016")),lines=list(col=c("darkred","blue"))),
col="darkred",
ylab=expression("mm","Wm"^-2,~degree~C,~degree~C,~degree~C,"ms"^-1,"ms"^-1,"ms"^-1),
ylim=list(c(0, 0.03),c(-0.1,-0.01),c(-0.1,-0.02),c(10,25),c(10,25),c(10,25),c(0,300),c(0,0.4)),
xlim=as.POSIXct(c("2015-03-01 00:00:00","2015-06-01 00:00:00")), more=TRUE)
bg <- trellis.par.get("background")
bg$col <- "transparent"
trellis.par.set("background", bg)
b=xyplot((stationseasonal[,c(1:8)]), layout=c(1,8),
scales=list(y=list(cex=0.7, rot=0)),
col="blue",
ylab=expression("mm","Wm"^-2,~degree~C,~degree~C,~degree~C,"ms"^-1,"ms"^-1,"ms"^-1),
xlim=as.POSIXct(c("2016-03-01 00:00:00","2016-06-01 00:00:00")))
a + as.layer(b, x.same=FALSE,axes = FALSE)
# WINTER
a=xyplot((stationseasonal[,c(1:8)]), layout=c(1,8), main="Winter",
scales=list(y=list(cex=0.7, rot=0)),
key=list(columns=2,text=list(lab=c("2014-2015","2015-2016")),lines=list(col=c("darkred","blue"))),
col="darkred",
ylab=expression("mm","Wm"^-2,~degree~C,~degree~C,~degree~C,"ms"^-1,"ms"^-1,"ms"^-1),
ylim=list(c(0, 0.04),c(-0.1,-0.02),c(-0.1,-0.05),c(5,20),c(5,20),c(5,20),c(0,200),c(0,0.4)),
xlim=as.POSIXct(c("2014-06-01 00:00:00","2014-09-01 00:00:00")), more=TRUE)
bg <- trellis.par.get("background")
bg$col <- "transparent"
trellis.par.set("background", bg)
b=xyplot((stationseasonal[,c(1:8)]), layout=c(1,8),
scales=list(y=list(cex=0.7, rot=0)),
col="blue",
ylab=expression("mm","Wm"^-2,~degree~C,~degree~C,~degree~C,"ms"^-1,"ms"^-1,"ms"^-1),
xlim=as.POSIXct(c("2015-06-01 00:00:00","2015-09-01 00:00:00")))
c=xyplot((stationseasonal[,c(1:8)]), layout=c(1,8),
scales=list(y=list(cex=0.7, rot=0)),
col="blue",
ylab=expression("mm","Wm"^-2,~degree~C,~degree~C,~degree~C,"ms"^-1,"ms"^-1,"ms"^-1),
xlim=as.POSIXct(c("2015-06-01 00:00:00","2015-09-01 00:00:00")))
a + as.layer(b, x.same=FALSE,axes = FALSE)
###################### TEST #########################
plot.new()
a=xyplot((stationplot[,c(1:7)]), layout=c(1,7), main="Spring",
scales=list(y=list(cex=0.7, rot=0)),
key=list(columns=2,text=list(lab=c("2014-2015","2015-2016")),lines=list(col=c("darkred","blue"))),
col="darkred",
ylab=expression("mm","Wm"^-2,~degree~C,~degree~C,"ms"^-1,"ms"^-1,"ms"^-1),
ylim=list(c(-0.02, 0.04),c(-0.08,0),c(-0.1,-0.02),c(10,30),c(10,30),c(10,30),c(0,400),c(0,0.4)),
xlim=as.POSIXct(c("2014-12-22 00:00:00","2015-04-01 00:00:00")), more=TRUE)
bg <- trellis.par.get("background")
bg$col <- "transparent"
trellis.par.set("background", bg)
b=xyplot((stationplot[,c(1:7)]), layout=c(1,7),
scales=list(y=list(cex=0.7, rot=0)),
col="blue",
ylab=expression("mm","Wm"^-2,~degree~C,~degree~C,"ms"^-1,"ms"^-1,"ms"^-1),
xlim=as.POSIXct(c("2015-12-22 00:00:00","2016-04-01 00:00:00")))
a + as.layer(b, x.same=FALSE,axes = FALSE)
## Auxiliary function to extract the year value of a POSIXct time
## index
Year <- function(x)format(x, "%Y")
xyplot(stationwindow, layout=c(1, ncol(stationwindow)), strip=FALSE,
scales=list(y=list(cex=0.8, rot=0)),
panel=function(x, y, ...){
## Values under the average highlighted with red regions
panel.xblocks(x, y<mean(y, na.rm=TRUE),
col = "indianred1",
height=unit(0.1, "npc"))
## Time series
panel.lines(x, y, col="royalblue4", lwd=0.5, ...)
## Label of each time series
panel.text(x[1], min(y, na.rm=TRUE),
names(stationwindow)[panel.number()],
cex=0.9, adj=c(0, 0), srt=90, ...)
## Triangles to point the maxima and minima
idxMax <- which.max(y)
panel.points(x[idxMax], y[idxMax],
col="black", fill="lightblue", pch=24)
idxMin <- which.min(y)
panel.points(x[idxMin], y[idxMin],
col="black", fill="lightblue", pch=25)
})
#Plot x acceleration (excluding wind direction)
dev.new()
xyplot(weatherzoo[,c(2,3,4,7)], layout=c(1, ncol(weatherzoo[,c(2,3,4,7)])), strip=FALSE,
scales=list(y=list(cex=0.8, rot=0)),
panel=function(x, y, ...){
## Time series
panel.lines(x, y, col="royalblue4", lwd=0.5, ...)
## Label of each time series
panel.text(x[1], min(y, na.rm=TRUE),
names(weatherzoo[,c(2,3,4,7)])[panel.number()],
cex=0.9, adj=c(0, 0), srt=90, ...)
## Triangles to point the maxima and minima
idxMax <- which.max(y)
panel.points(x[idxMax], y[idxMax],
col="black", fill="lightblue", pch=24)
idxMin <- which.min(y)
panel.points(x[idxMin], y[idxMin],
col="black", fill="lightblue", pch=25)
})
###########multiple plot
#Plot x,y,z acceleration with tides
dev.new()
plot(stationwindow[,c(1)])
par(new=TRUE)
plot(stationwindow[,c(4)], type="l", col="blue", xaxt="n", xlab="", ylab="")
axis(side=4)
mtext(side=4, line=4, 'fkfdsnkfv')
x <- rnorm(10)
y <- rnorm(10)
y2 <- rep(0, 10)
p1 <- xyplot(y ~ x, ylim=c, type="l")
p2 <- xyplot(y2 ~ x, pch=16, ylim=c(-5, 5), ylab=" ",type="l")
print(xyplot(station$X ~ station$Time), more=TRUE)
bg <- trellis.par.get("background")
bg$col <- "transparent"
trellis.par.set("background", bg)
print(p2)
print(c(a,b,layout=c(1,7)))
#check what this does
> library(mvtsplot)
> mvtsplot(stationwindow)
############################################################ January and july (24/10/16) ###################################################
#daily averages
#stationseasonal <- apply.daily((as.xts(stationplot)), colMeans)
#
stationplotraw <- xyplot(stationplot[,c(2,5,6,7,8)], layout=c(1,5), main="Raw station data",
scales=list(y=list(cex=0.7, rot=0)),
col="black", lwd=0.8,
xlim=as.POSIXct(c("2016-02-01 00:00:00","2016-03-01 00:00:00")),
ylim=list(c(-0.04,0.02),c(10,30),c(15,30),c(0,1000),c(0,6)),
ylab=expression("mm","Wm"^-2,~degree~C,~degree~C,"ms"^-2),more=TRUE)
update(stationplotraw, panel = function(...) {
panel.abline(h = 0, lty = "dotted", col = "red")
panel.xyplot(...)
})
# June july
a=xyplot((stationplot[,c(2,7)]), layout=c(1,2), main="Summer",
scales=list(y=list(cex=0.7, rot=0)),
key=list(columns=2,text=list(lab=c("2014-2015","2015-2016")),lines=list(col=c("darkred","blue"))),
col="dark grey",
ylab=expression("Wm"^-2,"ms"^-2),
ylim=list(c(-0.11,0.0),c(0,600)),
xlim=as.POSIXct(c("2016-06-01 00:00:00","2016-07-03 00:00:00")), more=TRUE)
bg <- trellis.par.get("background")
bg$col <- "transparent"
trellis.par.set("background", bg)
b=xyplot((stationplot[,c(3,8)]), layout=c(1,2),
scales=list(y=list(cex=0.7, rot=0)),
col="black",
ylab=expression("mm","ms"^-2),
ylim=list(c(-0.11,0.0),c(0,6)),
xlim=as.POSIXct(c("2016-06-01 00:00:00","2016-07-03 00:00:00")))
a + as.layer(b, x.same=TRUE, y.same=FALSE,axes = FALSE)
# march april
a=xyplot((stationplot[,c(2,7)]), layout=c(1,2), main="Summer",
scales=list(y=list(cex=0.7, rot=0)),
key=list(columns=2,text=list(lab=c("2014-2015","2015-2016")),lines=list(col=c("darkred","blue"))),
col="dark grey",
ylab=expression("Wm"^-2,"ms"^-2),
ylim=list(c(-0.11,0.05),c(0,800)),
xlim=as.POSIXct(c("2016-03-20 00:00:00","2016-04-20 00:00:00")), more=TRUE)
bg <- trellis.par.get("background")
bg$col <- "transparent"
trellis.par.set("background", bg)
b=xyplot((stationplot[,c(3,8)]), layout=c(1,2),
scales=list(y=list(cex=0.7, rot=0)),
col="black",
ylab=expression("mm","ms"^-2),
ylim=list(c(-0.11,0.05),c(0,6)),
xlim=as.POSIXct(c("2016-03-20 00:00:00","2016-04-20 00:00:00")))
a + as.layer(b, x.same=TRUE, y.same=FALSE,axes = FALSE)
# march april
a=xyplot((stationplot[,c(2,8)]), layout=c(1,2), main="Summer",
scales=list(y=list(cex=0.7, rot=0)),
key=list(columns=2,text=list(lab=c("2014-2015","2015-2016")),lines=list(col=c("darkred","blue"))),
col="dark grey",
ylab=expression("Wm"^-2,"ms"^-2),
ylim=list(c(-0.05,0.00),c(0,6)),
xlim=as.POSIXct(c("2016-03-20 00:00:00","2016-04-20 00:00:00")), more=TRUE)
bg <- trellis.par.get("background")
bg$col <- "transparent"
trellis.par.set("background", bg)
b=xyplot((stationplot[,c(3,6)]), layout=c(1,2),
scales=list(y=list(cex=0.7, rot=0)),
col="black",
ylab=expression("mm","ms"^-2),
ylim=list(c(-0.05,0.0),c(0,6)),
xlim=as.POSIXct(c("2016-03-20 00:00:00","2016-04-20 00:00:00")))
a + as.layer(b, x.same=TRUE, y.same=FALSE,axes = FALSE)
# march april
a=xyplot((stationplot[,c(2,5)]), layout=c(1,2), main="Summer",
scales=list(y=list(cex=0.7, rot=0)),
key=list(columns=2,text=list(lab=c("2014-2015","2015-2016")),lines=list(col=c("darkred","blue"))),
col="dark grey",
ylab=expression("Wm"^-2,"ms"^-2),
ylim=list(c(-0.11,0.05),c(10,22)),
xlim=as.POSIXct(c("2016-04-15 00:00:00","2016-04-16 00:00:00")), more=TRUE)
bg <- trellis.par.get("background")
bg$col <- "transparent"
trellis.par.set("background", bg)
b=xyplot((stationplot[,c(3,6)]), layout=c(1,2),
scales=list(y=list(cex=0.7, rot=0)),
col="black",
ylab=expression("mm","ms"^-2),
ylim=list(c(-0.11,0.05),c(10,22)),
xlim=as.POSIXct(c("2016-03-15 00:00:00","2016-04-16 00:00:00")))
a + as.layer(b, x.same=TRUE, y.same=TRUE,axes = FALSE)
# Define libraries required
library(zoo) 	# Provides access to function that forms a time series that is uneffected by missing/variable time data whislt also letting there be included time values for each day
library(EMD)	# Functions required for Empirical Mode Decomposition
library(hht)	# Hilbert-Huang Transform: Tools and Methods
library(Rssa)	# A Collection of Methods for Singular Spectrum Analysis
library(lattice) # Creates clean mulitvariate timeseries plots (zoo)
library(latticeExtra) # Added presentation options
library(TSA)	# Time series analysis package
# Close any open plots
graphics.off()
###IMPORTING DATA
# 10 minute weather station and accelerometer data 2014-2016
station <- read.csv("Omokoroa 10min data 2014 to Aug 2016.csv", header=T) # 10 min data
temp <- read.csv("internal temp.csv", header=T) # 10 min data
temp$hhnn<-as.POSIXct(paste(temp$dd.mm.yyyy, temp$hhnn), format="%d/%m/%Y %H:%M") # convert to time format "r" uses
temp <- temp[-c(1)]
colnames(temp)=c("Time","Internal Temperature")
station$Time<-as.POSIXct(paste(station$Date, station$Time), format="%d/%m/%Y %H:%M") # convert to time format "r" uses
station <- station[-c(1)]
colnames(station)=c("Time","X","Y","Z", "Air Temperature","Solar Radiation", "Soil Temperature","Wind Direction","Rainfall","Wind Speed") # changing coloums names
station <- merge(station,temp)
station <- station[,c(1,2,3,4,11,5,7,6,9,10,8)]
# convert mA acceleration into m/s2
station$X=((station$X*0.0000625)*9.81)
station$Y=((station$Y*0.0000625)*9.81)
station$Z=((station$Z*0.0000625)*9.81)
# Save or load binary version of converted data
save(station,file="stationbin.RData")
load("stationbin.RData")
###CONVERT TO REGULAR TIMESERIES (zoo)
# Convert one second data to zoo series
stationplot<-zoo(station[,2:11], station[,1])  #(Y[ROW,COLUMN], X[ROW, COLUMN])
stationplot2<-zoo(station[,2:4], station[,1])  #(Y[ROW,COLUMN], X[ROW, COLUMN])
summary(onesecondzoo)	# Display univariate statistics for columns in onesecondzoo (Index,X,Y,Z)
### EXTRACT SUBSETS OF DATA (defined intervals of data to reduce computation)
# Define a window of data and extract it
datawindow <- window(stationplot2, start = as.POSIXct("2016-07-03 00:00:00"), end = as.POSIXct("2016-07-11 00:00:00"))	# 10 minutes to match other data logger table
dataframewindow <- as.data.frame(datawindow)	# Converts time series into a dataframe that can be read by other programs such as Excel.
### PLOTTING RAW DATA
# basic plotting of various time series
title = paste(as.character(start(datawindow)),'to',as.character(end(datawindow)))
xyplot(datawindow, main=title, xlab="Time", ylab=expression('Acceleration m.s'^'-2')) # plot multivariate, separated time series plot
# Simple spectral analysis with 5% smoothing
dev.new()
spec <- spectrum(datawindow, spans=c(21,21,21), method="pgram")
spec$series <- "Accelerations"
plot(spec)
legend("right",legend=c("X axis","Y axis","Z axis"),lty=c(1,2,3),col=c(1,2,3))
# Simple Loess decomposition
dev.new()
plot(stl(ts(datawindow[,1], frequency=144), s.window="periodic"),main="X axis")
dev.new()
plot(stl(ts(datawindow[,2], frequency=144), s.window="periodic"),main="Y axis")
dev.new()
plot(stl(ts(datawindow[,3], frequency=144), s.window="periodic"),main="Z axis")
### EMD analysis
maximf <- 7	# Maximum number of empirical functions (IMFs) to fit (NULL forces the maximum possible)
# X axis
ttx <- seq(from=1,to=length(dataframewindow[,1]),by=1)	# Define a sequence of seconds
tryx <- emd(dataframewindow[,1],ttx,boundary="wave",max.imf=maximf)	# Undertake EMD
dev.new(width=6, height=10)		# Plot results in new window
par(mfrow=c(tryx$nimf+1, 1), mar=c(2,1,2,1))	# Define subplot regions for residual and number of IMFs found
rangeimfx <- range(tryx$imf)	# Define sclae of plot x asis to maximum range of IMF values
for(i in 1:tryx$nimf) {
plot(ttx, tryx$imf[,i], type="l", xlab="", ylab="", ylim=rangeimfx,
main=paste(i, "-th IMF X axis", sep="")); abline(h=0)
}		# Loop through to plot each IMF and then the residual
plot(ttx, tryx$residue, xlab="Seconds", ylab="", main="residual X axis", type="l", axes=FALSE); box()
# Y axis
tty <- seq(from=1,to=length(dataframewindow[,2]),by=1)
tryy <- emd(dataframewindow[,2],tty,boundary="wave",max.imf=maximf)
dev.new(width=6, height=10)
par(mfrow=c(tryy$nimf+1, 1), mar=c(2,1,2,1))
rangeimfy <- range(tryy$imf)
for(i in 1:tryy$nimf) {
plot(tty, tryy$imf[,i], type="l", xlab="", ylab="", ylim=rangeimfy,
main=paste(i, "-th IMF Y axis", sep="")); abline(h=0)
}
plot(tty, tryy$residue, xlab="Seconds", ylab="", main="residual Y Axis", type="l", axes=FALSE); box()
# Z axis
ttz <- seq(from=1,to=length(dataframewindow[,3]),by=1)
tryz <- emd(dataframewindow[,3],ttz,boundary="wave",max.imf=maximf)
dev.new(width=6, height=10)
par(mfrow=c(tryz$nimf+1, 1), mar=c(2,1,2,1))
rangeimfz <- range(tryz$imf)
for(i in 1:tryz$nimf) {
plot(ttz, tryz$imf[,i], type="l", xlab="", ylab="", ylim=rangeimfz,
main=paste(i, "-th IMF Z axis", sep="")); abline(h=0)
}
plot(ttz, tryz$residue, xlab="Seconds", ylab="", main="residual Z Axis", type="l", axes=FALSE); box()
### SSA analysis
# X axis
sxx <- ssa(dataframewindow[,1])		# Perform SSA on X axis to extract eigenvector orthogonal functions EOFs
summary(sxx)	# Display summary of annalysis
dev.new()
plot(sxx)		# Plot the proportion of variation accounted for by each EOF
dev.new()
plot(sxx,'vectors',idx=1:30)	# Plot the first 30 components based on EOFs
rxx <- reconstruct(sxx,groups=list(Trend=c(1),S1=c(2:13),S2=c(14:17),S3=c(18:30))) # Reconstruct into trend + 3 series (Need to tweak to select sensible EOFs)
dev.new()
plot(rxx, add.original = TRUE)  # Plot the reconstruction in a similar format to the EMD results
dev.new(width=6, height=10)
par(mfrow=c(4,1))
plot(rxx$S3,type="l",xlab="Seconds",ylab=expression('Acceleration m.s'^'-2'),main="S3 - X axis")
plot(rxx$S2,type="l",xlab="Seconds",ylab=expression('Acceleration m.s'^'-2'),main="S2 - X axis")
plot(rxx$S1,type="l",xlab="Seconds",ylab=expression('Acceleration m.s'^'-2'),main="S1 - X axis")
plot(rxx$Trend,type="l",xlab="Seconds",ylab=expression('Acceleration m.s'^'-2'),main="Trend - X axis")
# Y axis
syy <- ssa(dataframewindow[,2])
summary(syy)
dev.new()
plot(syy)
dev.new()
plot(syy,'vectors',idx=1:30)
ryy <- reconstruct(syy,groups=list(Trend=c(1),S1=c(2:15),S2=c(16:21),S3=c(22:30))) # Reconstruct into trend + 3 series
dev.new()
plot(ryy, add.original = TRUE)  # Plot the reconstruction
dev.new(width=6, height=10)
par(mfrow=c(4,1))
plot(ryy$S3,type="l",xlab="Seconds",ylab=expression('Acceleration m.s'^'-2'),main="S3 - Y axis")
plot(ryy$S2,type="l",xlab="Seconds",ylab=expression('Acceleration m.s'^'-2'),main="S2 - Y axis")
plot(ryy$S1,type="l",xlab="Seconds",ylab=expression('Acceleration m.s'^'-2'),main="S1 - Y axis")
plot(ryy$Trend,type="l",xlab="Seconds",ylab=expression('Acceleration m.s'^'-2'),main="Trend - Y axis")
# Z axis
szz <- ssa(dataframewindow[,3])
summary(szz)
dev.new()
plot(szz)
dev.new()
plot(szz,'vectors',idx=1:30)
rzz <- reconstruct(szz,groups=list(Trend=c(1),S1=c(2:9),S2=c(10:20),S3=c(21:30))) # Reconstruct into trend + 3 series
dev.new()
plot(rzz, add.original = TRUE)  # Plot the reconstruction
dev.new(width=6, height=10)
par(mfrow=c(4,1))
plot(rzz$S3,type="l",xlab="Seconds",ylab=expression('Acceleration m.s'^'-2'),main="S3 - Z axis")
plot(rzz$S2,type="l",xlab="Seconds",ylab=expression('Acceleration m.s'^'-2'),main="S2 - Z axis")
plot(rzz$S1,type="l",xlab="Seconds",ylab=expression('Acceleration m.s'^'-2'),main="S1 - Z axis")
plot(rzz$Trend,type="l",xlab="Seconds",ylab=expression('Acceleration m.s'^'-2'),main="Trend - Z axis")
test = rzz$S1+rzz$S2+rzz$S3
plot(test,type="l")
test2=ssa(test)
plot(test2)
plot(test2,'vectors',idx=1:30)
test3<-reconstruct(test2,groups=list(s1=c(1:4),s2=c(5:6),s3=c(7:26)))
plot(test3, add.original = FALSE)
dev.new(width=6, height=10)
par(mfrow=c(3,1))
plot(test3$s3,type="l",xlab="Seconds",ylab=expression('Acceleration m.s'^'-2'),main="S3 - Z axis")
plot(test3$s2,type="l",xlab="Seconds",ylab=expression('Acceleration m.s'^'-2'),main="S2 - Z axis")
plot(test3$s1,type="l",xlab="Seconds",ylab=expression('Acceleration m.s'^'-2'),main="S1 - Z axis")
### EMD analysis
maximf <- 7	# Maximum number of empirical functions (IMFs) to fit (NULL forces the maximum possible)
# X axis
ttx <- seq(from=1,to=length(dataframewindow[,1]),by=1)	# Define a sequence of seconds
tryx <- emd(dataframewindow[,1],ttx,boundary="wave",max.imf=maximf)	# Undertake EMD
dev.new(width=6, height=10)		# Plot results in new window
par(mfrow=c(tryx$nimf+1, 1), mar=c(2,1,2,1))	# Define subplot regions for residual and number of IMFs found
rangeimfx <- range(tryx$imf)	# Define sclae of plot x asis to maximum range of IMF values
for(i in 1:tryx$nimf) {
plot(ttx, tryx$imf[,i], type="l", xlab="", ylab="", ylim=rangeimfx,
main=paste(i, "-th IMF X axis", sep="")); abline(h=0)
}		# Loop through to plot each IMF and then the residual
plot(ttx, tryx$residue, xlab="Seconds", ylab="", main="residual X axis", type="l", axes=FALSE); box()
# Y axis
tty <- seq(from=1,to=length(dataframewindow[,2]),by=1)
tryy <- emd(dataframewindow[,2],tty,boundary="wave",max.imf=maximf)
dev.new(width=6, height=10)
par(mfrow=c(tryy$nimf+1, 1), mar=c(2,1,2,1))
rangeimfy <- range(tryy$imf)
for(i in 1:tryy$nimf) {
plot(tty, tryy$imf[,i], type="l", xlab="", ylab="", ylim=rangeimfy,
main=paste(i, "-th IMF Y axis", sep="")); abline(h=0)
}
plot(tty, tryy$residue, xlab="Seconds", ylab="", main="residual Y Axis", type="l", axes=FALSE); box()
# Z axis
ttz <- seq(from=1,to=length(dataframewindow[,3]),by=1)
tryz <- emd(dataframewindow[,3],ttz,boundary="wave",max.imf=maximf)
dev.new(width=6, height=10)
par(mfrow=c(tryz$nimf+1, 1), mar=c(2,1,2,1))
rangeimfz <- range(tryz$imf)
for(i in 1:tryz$nimf) {
plot(ttz, tryz$imf[,i], type="l", xlab="", ylab="", ylim=rangeimfz,
main=paste(i, "-th IMF Z axis", sep="")); abline(h=0)
}
plot(ttz, tryz$residue, xlab="Seconds", ylab="", main="residual Z Axis", type="l", axes=FALSE); box()
